%func:@使用解析模型计算圆柱磁铁的磁感应强度(B)
%inpt:@M:剩磁场
%     @(ao,l):磁铁参数(直径和高度)
%     @(xd,yd,zd)被计算的位置
%oupt:@B = [Bx;By;Bz]
%note:磁化方向变为了径向磁化
%date:2021-07-04 by xiyin,li in Hust;
%@Copyright 2021 Xiyin,li.
function B = PMCir3DByRadius(M,ao,l,position)
xd = position(1);
yd = position(2);
zd = position(3);
R = ao/2;
nu0 = 4*pi*10^(-7);
if(0)
    Bx1 = dblquad(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,0),-R,R,-l/2,l/2,1e-4)/(4*pi)*M;
    Bx2 = dblquad(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,1),-R,R,-l/2,l/2,1e-4)/(4*pi)*M;
    By1 = dblquad(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,2),-R,R,-l/2,l/2,1e-4)/(4*pi)*M;
    By2 = dblquad(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,3),-R,R,-l/2,l/2,1e-4)/(4*pi)*M;
    Bz1 = dblquad(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,4),-R,R,-l/2,l/2,1e-4)/(4*pi)*M;
    Bz2 = dblquad(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,5),-R,R,-l/2,l/2,1e-4)/(4*pi)*M;
else if(1)
        Bx1 = integral2(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,0),-R,R,-l/2,l/2,'RelTol',1e-4)/(4*pi)*M;
        Bx2 = integral2(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,1),-R,R,-l/2,l/2,'RelTol',1e-4)/(4*pi)*M;
        By1 = integral2(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,2),-R,R,-l/2,l/2,'RelTol',1e-4)/(4*pi)*M;
        By2 = integral2(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,3),-R,R,-l/2,l/2,'RelTol',1e-4)/(4*pi)*M;
        Bz1 = integral2(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,4),-R,R,-l/2,l/2,'RelTol',1e-4)/(4*pi)*M;
        Bz2 = integral2(@(x,z)PMCirFunc_R(R,x,z,xd,yd,zd,5),-R,R,-l/2,l/2,'RelTol',1e-4)/(4*pi)*M;
    end
end
B = [Bx1-Bx2;By1-By2;Bz1-Bz2];

%func:@使用解析模型计算球形磁铁时的被积函数
%inpt:@R:圆柱底面半径
%     @(x,z):被积平面的参数方程
%     @(xd,yd,zd)被计算的位置
%     @flag |0   |1   |2   |3   |4   |5
%     @axis |x↑ |x↓ |y↑ |y↓ |z↑ |z↓
%oupt:@difB = dif[Bx;By;Bz]
%note:
%date:2021-07-04 by xiyin,li in Hust;
%@Copyright 2021 Xiyin,li.
function dif_B = PMCirFunc_R(R,x,z,xd,yd,zd,flag)
%% 判断输入合法性
if(x>R)
    warnming("the x should NOT Greater than R")
end
%% 选择平面(右半柱面和左半柱面)
switch (mod(flag,2))
    case 0
        y = sqrt(R.^2-x.^2);
    case 1
        y = -sqrt(R.^2-x.^2);
end
R1 = sqrt((x-xd).^2+(y-yd).^2+(z-zd).^2);
%% 选择轴
switch (floor(flag/2)+1)
    case 1
        dif_B = (xd-x)./R1.^3;
    case 2
        dif_B = (yd-y)./R1.^3;
    case 3
        dif_B = (zd-z)./R1.^3;
    otherwise
        warning("The flag is NOT allowed!");
end
